Introducción a R como GIS

Vector

• Qué es un vector
• Cuantos tipos de vectores hay?
  • Puntos
  • Líneas
  • Polígonos
• Conocer la estructura interna de un vector
  • Número de vectores
  • Tabla de atributos

Raster

• qué es un raster?
• Qué es un extent?
• Cuáles son sus dimensiones?
• Saber qué es una proyección?
• Saber si un raster está proyectado o no?

Cómo leo un vector desde R.

El comando para leer un shapefile es readOGR; este se encuentra en la librería rgdal. Los argumentos del comando son:

• dsn: Es el folder donde tenemos nuestro archivo vector (shapefile)
• layer: El nombre de la capa que queremos leer

library(rgdal)

## rgdal: version: 1.3-6, (SVN revision 773)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
##  Linking to sp version: 1.3-1

library(sp)
map_dir <- "destdv1gw"
map_vec <-  readOGR(dsn=map_dir ,layer = "destdv1gw")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/Luis/Dropbox/CursoNichos2019/destdv1gw", layer: "destdv1gw"
## with 382 features
## It has 8 fields
## Integer64 fields read as strings:  COV_ COV_ID

Grafiquemos el mapa

plot(map_vec)

Explorando los atributos del vector

Uno de los comandos de más útilies para explorar los elementos que conforman a un objerto, es el comando str.

str(map_vec)

o bien para darnos solo una idea de su estrutura podemos escribir el nombre del objeto

map_vec

Los objetos espaciales en R son de la clase S4; para extraer los atributos de estos objetos (indexar) se utiliza el símbolo de **@**. Veamos como extraer la tabla de atributos del mapa map_vec

head(map_vec@data)

Número de vectores o poligonos

length(map_vec@plotOrder)

## [1] 382

Extraer un polígono del shapefile

puebla <- map_vec[map_vec@data$ENTIDAD=="PUEBLA",]
puebla <-  map_vec[map_vec@data$ENTIDAD %in% "PUEBLA",]
plot(puebla)

Guardar un shapefile

Ahora guardamos como un shapefile separado, el polígono

writeOGR(obj = puebla,layer = "puebla",
         dsn = "puebla",driver = "ESRI Shapefile")

Filtrar los puntos que caen dentro de un polígono

Generemos unos puntos aleatorios definidos en el extent del mapa; el comando nombre_objeto@bbox nos muestra el extent de nuestro shapefile.

map_vec@bbox

##          min       max
## x -118.36603 -86.71074
## y   14.53401  32.71877

# Coordenadas aleatorias
rand_x <- runif(n =300, map_vec@bbox[1,1],max = map_vec@bbox[1,2])
rand_y <- runif(n =300, map_vec@bbox[2,1],max = map_vec@bbox[2,2])
coords_rand <- data.frame(longitude=rand_x,latitude=rand_y)

# Cual es el CRS de map_vec
crs_mexico <- map_vec@proj4string@projargs
sp_df_rand <- SpatialPoints(coords = coords_rand,proj4string =CRS(crs_mexico))

plot(map_vec)
plot(sp_df_rand,add=T)

data_poly_all <- sp::over( sp_df_rand , map_vec , fn = NULL)
en_poligono_index <- which(!is.na(data_poly_all[,1]))
p_en_poligono <- sp_df_rand[en_poligono_index ,]
plot(map_vec)
plot(p_en_poligono,add=T,col="blue")

data_en_poli_extract <- data_poly_all[en_poligono_index,]
coords_en_poli <- sp_df_rand@coords[en_poligono_index,]
data_en_poli_extract_df <- cbind(coords_en_poli,data_en_poli_extract)
sp_DF_en_poli_extract <- SpatialPointsDataFrame(data_en_poli_extract_df[,c("longitude","latitude")],
                                                data_en_poli_extract_df)
#writeOGR(sp_DF_en_poli_extract,
#         dsn = map_dir,layer = "puntos_aleatorios_mexico",
#         driver = "ESRI Shapefile")

knitr::kable(head(sp_DF_en_poli_extract@data))

#plot(map_vec)
#cols <- unique(as.numeric(sp_DF_en_poli_extract@data$ENTIDAD))
#color_fun <- function(x) return()
#plot(p_en_poligono,add=T,col=sp_DF_en_poli_extract@data$ENTIDAD)

Reproyectar un polígono

Con R podemos hacer reproyección de coordenadas, podemos pasar de coordenadas geográficas (logitude, latitud) a coordenadas utm (metrosal nivel del mar).

Veamo como proyectar nuesto mapa de Puebla a coordenadas utm

puebla_reproj <- spTransform(puebla,CRS("+proj=utm +zone=10+datum=WGS84"))
plot(puebla_reproj,axes=TRUE)

Regresemos a coordenadas geográficas

puebla_reproj_geo <- spTransform(puebla_reproj,CRS("+proj=longlat +ellps=WGS84"))
plot(puebla_reproj_geo,axes=TRUE)

Límites territoriales

• GADM: Datos de límites territoriales de países

Datos climáticos

• worldclim: Datos climáticos

Datos Digitales de Elevación

• SRTM: Datos Digitales de Elevación SRTM 90m.
• mapzen: Datos Digitales de Elevación aprox ~ 9m.

library(raster)
wc10 <- raster::getData('worldclim', var='bio', res=10)
plot(wc10)

Cómo leer un raster

El comando para leer 1 solo raster es raster

library(raster)
bio2 <- raster("wc10/bio2.bil")
plot(bio2)

Cómo leer un conjunto (stack) de capas

A veces (casi siempre) será necesario leer un conjunto de capas que tienen la misma resolución y extent (i.e. WorldClim layers). El comando stack de raster permite crear stack de capas; para ello es necesario indicar con en un vector los path o rutas a la carpeta que contiene nuestras capas.

Veamos como leer o almacenar en la memoria de R el conjunto de capas de WorldClim layers

Primero en guardamos en una variable los path a nustras capas

paths_capas <- list.files("wc10/",
                          pattern = "*.bil$",full.names = TRUE)

Usamos el comando stack

bios_wc <- raster::stack(paths_capas)

Graficamos las primer 5 capas

plot(bios_wc[[1:5]])

notemos que utilice la misma sintaxis que utilizamos para indexa listas…

Hacer un crop

La función crop utiliza el extent de un shapfile para hacer recortar un raster o un stack. Veamos un recorte de las capas de WorldClim con el poligono de puebla.

puebla_ras <- raster::crop(x = bios_wc,y = puebla)
plot(puebla_ras[[1]])
plot(puebla,add=T)

Hacer una mascara

Al igual que crop, mask hace un recorte pero conserva la forma del poligono.

puebla_ras_forma <- raster::mask(puebla_ras, puebla)
plot(puebla_ras_forma,add=T)

Guardar un raster

El comando para gurardar un raster es wirteRaster y está en el paquete raster

# Guardar el crop de puebla
writeRaster(puebla_ras[[1]],"puebla_bio2.tif",overwrite=TRUE)

Guardar un stack completo

Mostraremos como guardar un stack completo en una carpeta

# Creamos la carpeta
dir.create("Puebla_bios")
lapply(names(puebla_ras_forma), function(x){
  writeRaster(puebla_ras_forma[[x]], paste0("Puebla_bios/",
                        x,".asc"),overwrite=TRUE)})

Convertir a puntos un stack

stack_df <- data.frame(rasterToPoints(puebla_ras_forma))
head(stack_df)

Extraer valores de un raster o un stack

Podemos extraer los valores de un stack utilizando coordenadas.

plot(bios_wc[[1]])
points(coords_rand,pch=".",col="red")

coord_rand_extract <- data.frame(extract(bios_wc,coords_rand))
coord_rand_extract <- na.omit(coord_rand_extract)
knitr::kable(head(coord_rand_extract))

Tabla de correlaciones…

knitr::kable(round(cor(coord_rand_extract),2))

Reproyectar un raster

Al igual que con los vectores, podemos reproyectar a los rasters. Lo anterior se puede hacer mediante el comando projectRaster de raster

puebla_bio1 <- puebla_ras[[1]]
puebla_bio1_reproj <- projectRaster(puebla_bio1,
                                    crs="+proj=utm +zone=10+datum=WGS84")
plot(puebla_bio1_reproj)
plot(puebla_reproj,add=T)

Regresemos a geográficas

puebla_bio1_geo <- projectRaster(puebla_bio1, 
                                 crs="+proj=longlat +ellps=WGS84")
plot(puebla_bio1_geo)
plot(puebla,add=T)

Descarga de datos de presencia spocc

library(spocc)

jaguar_occs <- occ("Panthera onca",
                    from = "gbif",
                    limit = 1500,
                    has_coords = T)
jaguar_occsDF <- occ2df(jaguar_occs)
map %>% leaflet::addCircleMarkers(lng =  jaguar_occsDF$longitude,
                                  lat = jaguar_occsDF$latitude)

Un poco más profesional

# install.packages("wicket")
library(wicket)
mex <- wrld_simpl[wrld_simpl$NAME %in% c("Mexico"),]
raster::extent(mex)

## class       : Extent 
## xmin        : -118.4042 
## xmax        : -86.7014 
## ymin        : 14.55055 
## ymax        : 32.71846

jaguar_mexExtent <- occ(query = 'Panthera onca',
                        geometry =c(-118.40417,14.55055,
                                    -86.70140,32.71846))
jaguar_mexExtentDF <- occ2df(jaguar_mexExtent)
map %>% leaflet::addCircleMarkers(lng =  jaguar_mexExtentDF$longitude,
                                  lat = jaguar_mexExtentDF$latitude)

Aún más pro

jaguar_mex <- occ(query = 'Panthera onca', 
                  geometry = wicket::sp_convert(mex))
jaguar_mexDF <- occ2df(jaguar_mex)
map %>% leaflet::addCircleMarkers(lng =  jaguar_mexDF$longitude,
                                  lat = jaguar_mexDF$latitude)

Revisen su directorio de trabajo